In-Class Exercise 5

A short description of the post.

Installing and Loading the R package

• maptools is for creating spatial objects
• sf is for importing geospatial data
• raster is for creating raster object
• spatstat is for performing SPPA
• tmap is for plotting thematic maps
• tidyverse is for data handling

packages = c('maptools', 'sf', 'raster', 'spatstat', 'tmap', 'tidyverse', 'plotly', 'ggthemes')
for (p in packages){
if(!require(p, character.only = T)){
install.packages(p)
}
library(p,character.only = T)
}

Importing the Geospatial Data

Importing Shapefile using st_read() of sf package.The output object is in tibble sf object class.

mpsz_sf <- st_read(dsn = "data/shapefile", layer = "MP14_SUBZONE_WEB_PL")

Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `C:\darrylkwok\IS415_blog\_posts\2021-09-13-in-class-exercise-5\data\shapefile' 
  using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21

mpsz_sf

Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21
First 10 features:
   OBJECTID SUBZONE_NO       SUBZONE_N SUBZONE_C CA_IND
1         1          1    MARINA SOUTH    MSSZ01      Y
2         2          1    PEARL'S HILL    OTSZ01      Y
3         3          3       BOAT QUAY    SRSZ03      Y
4         4          8  HENDERSON HILL    BMSZ08      N
5         5          3         REDHILL    BMSZ03      N
6         6          7  ALEXANDRA HILL    BMSZ07      N
7         7          9   BUKIT HO SWEE    BMSZ09      N
8         8          2     CLARKE QUAY    SRSZ02      Y
9         9         13 PASIR PANJANG 1    QTSZ13      N
10       10          7       QUEENSWAY    QTSZ07      N
        PLN_AREA_N PLN_AREA_C       REGION_N REGION_C
1     MARINA SOUTH         MS CENTRAL REGION       CR
2           OUTRAM         OT CENTRAL REGION       CR
3  SINGAPORE RIVER         SR CENTRAL REGION       CR
4      BUKIT MERAH         BM CENTRAL REGION       CR
5      BUKIT MERAH         BM CENTRAL REGION       CR
6      BUKIT MERAH         BM CENTRAL REGION       CR
7      BUKIT MERAH         BM CENTRAL REGION       CR
8  SINGAPORE RIVER         SR CENTRAL REGION       CR
9       QUEENSTOWN         QT CENTRAL REGION       CR
10      QUEENSTOWN         QT CENTRAL REGION       CR
            INC_CRC FMEL_UPD_D   X_ADDR   Y_ADDR SHAPE_Leng
1  5ED7EB253F99252E 2014-12-05 31595.84 29220.19   5267.381
2  8C7149B9EB32EEFC 2014-12-05 28679.06 29782.05   3506.107
3  C35FEFF02B13E0E5 2014-12-05 29654.96 29974.66   1740.926
4  3775D82C5DDBEFBD 2014-12-05 26782.83 29933.77   3313.625
5  85D9ABEF0A40678F 2014-12-05 26201.96 30005.70   2825.594
6  9D286521EF5E3B59 2014-12-05 25358.82 29991.38   4428.913
7  7839A8577144EFE2 2014-12-05 27680.06 30230.86   3275.312
8  48661DC0FBA09F7A 2014-12-05 29253.21 30222.86   2208.619
9  1F721290C421BFAB 2014-12-05 22077.34 29893.78   6571.323
10 3580D2AFFBEE914C 2014-12-05 24168.31 30104.18   3454.239
   SHAPE_Area                       geometry
1   1630379.3 MULTIPOLYGON (((31495.56 30...
2    559816.2 MULTIPOLYGON (((29092.28 30...
3    160807.5 MULTIPOLYGON (((29932.33 29...
4    595428.9 MULTIPOLYGON (((27131.28 30...
5    387429.4 MULTIPOLYGON (((26451.03 30...
6   1030378.8 MULTIPOLYGON (((25899.7 297...
7    551732.0 MULTIPOLYGON (((27746.95 30...
8    290184.7 MULTIPOLYGON (((29351.26 29...
9   1084792.3 MULTIPOLYGON (((20996.49 30...
10   631644.3 MULTIPOLYGON (((24472.11 29...

Projection is in SVY21.

Importing the aspatial data

read_rds() of readr package is used instead of readRDS() of base R. This is because the output of read_rds() is a tibble object.

CHAS <- read_rds("data/rds/CHAS.rds")
childcare <- read_rds("data/rds/childcare.rds")

Note that there are some issues found in the childcare dataframe, because the Lat and Lng attributes should be in numeric data type. The coordinate fields seems to be in decimal degree. Therefore, this results in having assumed the WGS84 referencing system.

st_crs(CHAS)

Coordinate Reference System: NA

Converting the aspatial data into sf objects

CHAS_sf <- st_as_sf(CHAS, coords = c("X_COORDINATE", "Y_COORDINATE"), crs = 3414)

Note: st_as_sf accepts coordinates that are in character data type, and converts it into the appropriate data type.

We see that it is in decimal degree in the childcare data, therefore we know that it is 4326 projection system and we will transfrom it to 3414

childcare_sf <- st_as_sf(childcare, coords = c("Lng", "Lat"), crs = 4326) %>%
  st_transform(crs=3414)

Plotting for initial review of the Geospatial data

This is how we can check the projection system. If the points are plotted on somewhere random, it means that the projection system is wrong.

All the darker points show that you have more than 1 childcare centers

tmap_mode("view")

tm_shape(childcare_sf) +
  tm_dots(alpha = 0.4, col = "blue", size = 0.05) +
tm_shape(CHAS_sf) + 
  tm_dots(alpha = 0.4, col = "red", size = 0.05)

Geospatial Data Wrangling

Convert sf objects to Spatial Classes

childcare <- as_Spatial(childcare_sf)
CHAS <- as_Spatial(CHAS_sf)
mpsz <- as_Spatial(mpsz_sf)

Convert Spatial Data into Spatial Objects

Using *as.ppp() of maptools** package

childcare_sp <- as(childcare, "SpatialPoints")
CHAS_sp <- as(CHAS, "SpatialPoints")
mpsz_sp <- as(mpsz, "SpatialPolygons")

Convert Spatial Objects into ppp objects

Using as.ppp() of maptools package

childcare_ppp <- as(childcare_sp, "ppp")
CHAS_ppp <- as(CHAS_sp, "ppp")

Removing duplicate points using jitter

If you check childcare_ppp and CHAS_ppp for duplicates using the any() and duplicated() function, you will see that it will return TRUE

childcare_ppp_jit <- rjitter(childcare_ppp, retry=TRUE, nsim=1, drop=TRUE)

any(duplicated(childcare_ppp_jit))

[1] FALSE

CHAS_ppp_jit <- rjitter(CHAS_ppp, retry=TRUE, nsim=1, drop=TRUE)

any(duplicated(CHAS_ppp_jit))

[1] FALSE

Extracting Punggol Planning Area

Need a comma at the end to select all the columns since the statement before the comma is to select the rows

pg <- mpsz[mpsz@data$PLN_AREA_N == "PUNGGOL",]

Convert pg, which is a SpatialPolygonDataFrame into SpatialPolygonObject, called pg_sp*

pg_sp <- as(pg, "SpatialPolygons")

Convert pg_sp into a owin object

pg_owin <- as(pg_sp, "owin")

Extract spatial points within owin object

Extract out the childcares and clinics that are within punggol.

If you plot childcare_ppp_jit, you will only get the points. After extracting the points using the owin object, you will get the Punggol’s Boundary and the points within Punggol

childcare_pg <- childcare_ppp_jit[pg_owin]
CHAS_pg <- CHAS_ppp_jit[pg_owin]

plot(childcare_pg)

L-function

L_childcare <- envelope(childcare_pg, Lest, nsim=99, rank=1, global=TRUE)

Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,
71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.

Done.

Code chunk for plotting interactive L-function

title <- "Pairwise Distance: L function"

Lcsr_df <- as.data.frame(L_childcare)

colour=c("#0D657D","#ee770d","#D3D3D3")
csr_plot <- ggplot(Lcsr_df, aes(r, obs-r))+
  # plot observed value
  geom_line(colour=c("#4d4d4d"))+
  geom_line(aes(r,theo-r), colour="red", linetype = "dashed")+
  # plot simulation envelopes
  geom_ribbon(aes(ymin=lo-r,ymax=hi-r),alpha=0.1, colour=c("#91bfdb")) +
  xlab("Distance r (m)") +
  ylab("L(r)-r") +
  geom_rug(data=Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,], sides="b", colour=colour[1])  +
  geom_rug(data=Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,], sides="b", colour=colour[2]) +
  geom_rug(data=Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,], sides="b", color=colour[3]) +
  theme_tufte()+
  ggtitle(title)

text1<-"Significant clustering"
text2<-"Significant segregation"
text3<-"Not significant clustering/segregation"

# the below conditional statement is required to ensure that the labels (text1/2/3) are assigned to the correct traces
if (nrow(Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,])==0){ 
  if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){ 
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text3, traces = 4) %>%
      rangeslider() 
  }else if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){ 
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text2, traces = 4) %>%
      rangeslider() 
  }else {
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text2, traces = 4) %>%
      style(text = text3, traces = 5) %>%
      rangeslider() 
  }
} else if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){
  if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text1, traces = 4) %>%
      rangeslider() 
  } else{
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text1, traces = 4) %>%
      style(text = text3, traces = 5) %>%
      rangeslider()
  }
} else{
  ggplotly(csr_plot, dynamicTicks=T) %>%
    style(text = text1, traces = 4) %>%
    style(text = text2, traces = 5) %>%
    style(text = text3, traces = 6) %>%
    rangeslider()
  }